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Abstract. We study stellar convection using a local three-dimensional MHD model, with which we investigate 
the influence of rotation and large-scale magnetic flelds on the turbulent momentum and heat transport and 
their role in generating large-scale flows in stellar convection zones. The former is studied by computing the 
turbulent velocity correlations, known as Reynolds stresses, the latter by calculating the correlation of velocity 
and temperature fluctuations, both as functions of rotation and latitude. We find that the horizontal correlation, 
Qsti), capable of generating horizontal differential rotation, attains significant values and is mostly negative in 
the southern hemisphere for Coriolis numbers exceeding unity, corresponding to equatorward fiux of angular 
momentum. This result is also in accordance with solar observations. The radial component Qr4, is negative for 
slow and intermediate rotation indicating inward transport of angular momentum, while for rapid rotation, the 
transport occurs outwards. Parametrisation in terms of the mean-field A-effect shows qualitative agreement with 
the turbulence model of Kichatinov & Riidiger (1993 ) for the horizontal part H cx , whereas for the vertical 
A-effect, V OC ^^^^ , cL greement only for intermediate rotation exists. The A-coefficients become suppressed in the 
limit of rapid rotation, this rotational quenching being stronger and occurring with slower rotation for the V 
component than for H. We have also studied the behaviour of the Reynolds stresses under the influence of a 
large-scale azimuthal magnetic held of varying strength. We flnd that the stresses are enhanced by the presence 
of the magnetic held for fleld strengths up to and above the equipartition value, without signiflcant quenching. 
Concerning the turbulent heat transport, our calculations show that the transport in the radial direction is most 
efficient at the equatorial regions, obtains a minimum at midlatitudes, and shows a slight increase towards the 
poles. The latitudinal heat transport does not show a systematic trend as a function of latitude or rotation. 
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1. Introduction 

The Sun is rotating differentially; its internal rotation pro- 
file is known from helioseismic inversions (e.g. Schou et al. 
119981 Thompson et al. I2003|l . There also exists a merid- 
ional large-scale flow, although the magnitude and struc- 
ture of this meridional circulation pattern is rather poorly 
known observationally (e.g. Stix 1200^ . For more active 
rapidly rotating stars, photometric (e.g. Hall |199T| Henry 
et al. TO9 5 Donahue et al. 11996(1 and spectroscopic (e.g. 
Collier Cameron et al. IMT^ Reiners & Schmitt I2003al 
I2003b|l observations have been used to deduce the surface 
difi^erential rotation. These observations show that the rel- 
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ative differential rotation A^l/^l strongly decreases as a 
function of rotation, whilst the absolute one, Afi, remains 
roughly constant. 

The solar (and to some extent stellar) rotation has 
been studied using three main approaches: local and global 
convection models and mean-field calculations. The global 
convection models, which have been in use for more than 
two decades (e.g. Gilman I1977f) improving in the cover- 
age of spatial scales, have only very recently started to 
show reasonable agreement with the helioseismic observa- 
tions (e.g. Robinson & Chan lWTl Brun & Toomre OUU^ . 
The clear advantage of these models is the fact that the 
Reynolds stresses, meridional flows, differential rotation 
and their non-linear interactions can be studied from the 
same model. These calculations, however, are computa- 
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tionally expensive, making a comprehensive study of dif- 
ferent types of stars with varying rotation and depths 
of convection zones difficult to undertake. Furthermore, 
the spatial scales covered are not yet adequately small to 
consistently model the small-scale turbulence, which may 
have an important influence on the dynamics of the large- 
scale flow. 

A theory of the role of anisotropic turbulence on an- 
gular momentum transport and large-scale motions was 
formulated by Riidiger (1 977l see also liiidiger [19891 here- 
after R89, for a thorough discussion), in the form of the 
so-called A-effect: in a rotating anisotropic fluid additional 
Reynolds stresses arise, which are proportional to H. itself, 
therefore non-vanishing even for rigid rotation; in stel- 
lar convection zones the radial density stratification pro- 
vides a natural source of the required anisotropy. These 
turbulent stresses, as an ensemble average, lead to an- 
gular momentum transport both horizontally and verti- 
cally, giving rise to non-uniform rotation profiles. This 
can be quantified, under the mean-field approximation, 
as the components of the A-tensor in the series expansion 
Qij = A^jk^k +-N'i]kid (uk) /dxi. However, a major diffi- 
culty remains: how to determine the magnitude and form 
of the tensor components for the convectively turbulent 
flow, for which no complete theory exists. 

The small-scale turbulent motions can also take part in 
the generation of meridional flows, due to the rotationally 
induced latitude dependence of the turbulent heat flux, 
measured by the correlation tensor {u[T') (e.g. R89). The 
meridional component of the Reynolds stress tensor, Qre^ 
can also have an influence on the meridional flows (e.g. 
R89) , although the significance of this component has not 
yet been fully explored. Any meridional flow generated 
will take part in the angular momentum transport, ca- 
pable of driving horizontal differential rotation; therefore, 
its investigation is of importance in determining the stellar 
rotation profiles. 

In most cases the parametrisations for the turbulent 
quantities are derived under some simplifying assump- 
tions, most notably the first order smoothing approxi- 
mation (FOSA), also called the second order correlation 
approximation (SOCA) or the quasi-linear approach, in 
which no terms higher than second order in the linearised 
equations for the fiuctuating fields are taken into account. 
This is valid in the hydrodynamic case if the Strouhal 
number, which is the ratio of the coherence time to the 
turnover time, is small. This condition is not very well 
satisfied in the stellar environment (e.g. Stix 2002, sec 
also Brandenburg et al. .2004 ). Nevertheless, models ap- 
plying FOSA have been used to derive the A-tensor coeffi- 
cients (recently e.g. by Kichatinov & Riidiger 1X9981 here- 
after KR93), the diffusive part of the turbulent stresses 
(Kitchatinov et al. I1994a|l . and to investigate the be- 
haviour of the A-coefficients in the magnetised regime 
(Kitchatinov et al. .1994b ). 

The resulting formulations of the transport coefficients 
have been used to study stellar rotation with mean-field 
models, in which approach, in the simplest form, only the 



mean-field momentum (Reynolds) equation is solved for 
(e.g. Kitchatinov et al. ll994aj) . More sophisticated models 
(e.g. Brandenburg et al. ll992|l work with the MHD equa- 
tions taking into account the thermodynamics as well. The 
advantage of these models is the low computational cost, 
allowing for the exploration of large parameter regimes. 
Using this approach, rotation laws of the young solar- 
like stars and the Sun, as well as some other stars (e.g. 
Kichatinov & Riidiger EMI HMSIIMI Kiiker et al.EMSI 
Kiiker & Stix .2001i) have been calculated. Along a different 
line of study, Rieutord et al. H1994|) demonstrated that the 
differential rotation from a direct convection calculation in 
spherical coordinates agrees with a mean-field model using 
a A-effect derived from direct calculation. 

Another way to derive the turbulent transport coeffi- 
cients, bypassing the difficulties of solving the full prob- 
lem in a spherical domain and the probably oversimplified 
analytical techniques, is to restrict the calculations to a 
small Cartesian domain, located at different latitudes and 
depths in the stellar convection zone. Following this kind 
of an approach, Hathaway & Somerville l|1988|l . Pulkkinen 
et al. 14993,) . Brummell et al. H1998|l . and, in a more com- 
prehensive study, Chan H2001() . have calculated the off- 
diagonal components of the Reynolds stress tensor. In the 
study of Pulkkinen et al. (jl993|l the turbulent heat trans- 
port was investigated as well. 

The investigation of solar and stellar rotation plays 
an important role in understanding the stellar magnetic 
activity as well. This is because the solar and stellar mag- 
netic fields are believed to arise due to the inductive ac- 
tion of turbulent motions and differential rotation, which 
is generally called the turbulent dynamo mechanism. In 
the special case of the Sun, for which the rotation profile 
is an observed quantity, the dynamo models succesfuUy 
reproduce the solar magnetic activity patterns only if a 
relatively strong meridional circulation is included (e.g. 
Bonnano et al. I2002|l : since the observational knowledge 
of such motions is still rather poor, it is important to 
study them from numerical models both on global and 
local scales. For the more active rapidly rotating stars, 
observational information is limited to the deduction of 
the surface differential rotation, giving an even more im- 
portant role for the numerical investigations. 

In this study, we investigate the characteristics and 
correlations in a turbulent convective flow, using a local 
three-dimensional MHD model of Caunt & Korpi (2001), 
on top of which we have implemented a three-layer setup 
consisting of an overshoot layer on the bottom, a convec- 
tively unstable layer in the middle, and a narrow cooling 
layer on top of the domain. As in the previous studies, 
we calculate the off-diagonal components of the Reynolds 
stresses and radial and meridional components of the tur- 
bulent heat transport tensor, but extend our calculations 
further to the rapid rotation end (up to Coriolis number 
15), and also investigate the infiuence of imposed magnetic 
fields on the turbulent dynamics. As our model includes 
a stable overshoot layer below the convectively unstable 
one, we also report on the relevance of this layer for the 
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turbulent dynamics, and in addition on the dependence of 
the overshooting efficiency as a function of rotation and 
latitude. 

The remainder of the paper is organised as follows: in 
Sect. 12 the model is described, in Sect.Othe parameters 
of the calculations are given and the methods used are 
discussed, and in Sect.0] the mean-field theory of the A- 
effect and turbulent heat transport is introduced. Sects. 
and^lpresent the results and the conclusions, respectively. 

2. The model 

The computational domain in our calculations is a rectan- 
gular box, situated on the southern hemisphere of a star 
at a latitude 6, see Fig. ^ The coordinates are chosen 
so that X, y, and z correspond to the south-north, west- 
east, and radially inward directions, respectively. When 
analysing the results it is important to keep in mind 
that our {x,y,z) corresponds to {—9,(p,—r) in spherical 
polar coordinates, where 6 is the colatitude. The angu- 
lar velocity as function of latitude can now be written 
as n = r2(cos0e3; — sinOcz). We solve a set of MHD- 
equations 



dA 

'dt 
dhip 

dt 
du 

'at 



de 

'Ft 



= u X B 'i]poJ , (1) 

= -(tt- V)lnp + V-M, (2) 

= -(m • V)m - -Vp -2n X u + g + 
P 

+ -JxB + -V-a, (3) 
P P 

V 1 
= -{u ■ V)e - -(V • It) + -V • (xpVe) + 

P P 

"i^Tvisc ~t- Tjoulc rcool i (4) 



where A is the vector potential, u the velocity, B = V x A 
the magnetic field, J = V x B/ the current density, p 
the mass density, p the pressure, g = ge^ the (constant) 
gravity, and e = cyT the internal energy per unit mass, rj 
and v are the magnetic diffusivity and kinematic viscosity, 
respectively, /iq the vacuum permeability, and x the ther- 
mal diffusivity. The equation of state is that of an ideal 
gas 

p = pe(7-l), (5) 

where the ratio of the specific heats 7 = cp/cy = 5/3. 
cr — 2pvS is the stress tensor where 



1 / dui 9m, 



2 V dxj dxi 



- SijV ■ u 



(6) 



The terms responsible for viscous and Joule heating can 
be written as 



Tvisc 2uSij 



dxi 



-L Joule — *J 

p 



(7) 
(8) 



We use a narrow cooling layer on top of the convection 
zone, cooled with a term 



1^000! 



1 



cool 



■/(z)(e- eo) , 



(9) 



where icooi is a cooling time, chosen to be short enough 
for the upper boundary to stay isothermal, f{z) a func- 
tion which vanishes everywhere else but in the interval 

< z < zi, and cq — e(zo) the value of internal energy 
at the top of the box. This parametrisation mimics the 
radiative losses at the stellar surface and, although still 
rather simple in comparison to the real surface layers of 
stars, works as a more realistic boundary condition and 
stabilises the numerics better than just imposing a con- 
stant temperature at the boundary. 

We adopt periodic boundary conditions in the horizon- 
tal directions, and closed stress free boundaries at the top 
and bottom. The temperature is kept fixed at the top of 
the box and a constant heat flux is applied at the bottom 



dz 



duy 
dz 



e{zo) 
de 
dz 





eo , 



(7-1)(to3 + 1) ' 



at z = zq, Z3 ; (10) 

(11) 
(12) 



where is the polytropic index of the lower overshoot 
layer. 

We adopt the same dimensionless quantities and box 

mm 



dimensions as Brandenburg et al 
et al. PMIIM?^ . Thus by setting 

d = Pq = g = po = cp = 1 , 



and Ossendrijver 



length is measured with respect to the depth of the un- 
stable layer, d — Z2 — zi, and density in units of the initial 
value at the bottom of the convectively unstable layer, pq. 
Time is measured in units of the free fall time, ^/djg, ve- 
locity in units of ^/dg, magnetic field in terms of ^dpopog, 
and entropy in terms of cp. The box has horizontal di- 
mensions Lx = Ly = 4, and Lz = 2 in the vertical direc- 
tion. The vertical coordinates are set to be (zq, zi, Z2, 23) 
— (—0.15, 0, 1, 1.85), where zq and Z3 are the top and bot- 
tom of the box, and zi and Z2 are the boundaries between 
stable and unstable layers. 

The dimensionless parameters controlling the calcula- 
tions are the kinematic and magnetic Prandtl numbers Pr 
and Pm, and the Taylor, Rayleigh, and Chandrasekhar 
numbers, denoted by Ta, Ra, and Ch, respectively. The 
relative importances of thermal and magnetic diffusion 
against the kinetic one are measured by the kinetic and 
magnetic Prandtl numbers 



Pr = 



Pm 



Xo 

V 

77 



(13) 
(14) 



where xo is the reference value of the thermal diffusivity, 
taken from the middle of the unstably stratified layer. In 
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Fig. 1. Sketch of the model setup. Left: the coordmate system is set up so that x-axis points from south to north, 
y-axis from west to east, and z-axis radially inwards. The angle O is the latitude, i.e. the angular distance of the centre 
of the box from the equator. Right: the domain is divided into three parts, an upper cooling layer, the convectively 
unstable layer, and the lower overshoot region. The coordinates zq and Z3 denote the upper and lower boundaries of 
the box, whereas zi and Z2 are the boundaries between stable and unstable layers. 



the present study, v and ry are constants. We set Pr — 
0.4 which is a compromise between computational time 
(large diffusivity enforces a smaller time step) and realistic 
physics (in the Sun Pr <C 1). In the magnetic runs we set 
Pm = 1 because of a similar compromise between realistic 
physics and the need to keep Rm reasonably large. 
Rotation is measured by the Taylor number 



Ta 



(15) 



A related dimensionless quantity is the Coriolis number, 
which is twice the inverse of the Rossby number, Co = 
2Qt, where r = d/ut is the convective turnover time, 
and ut — {u''^)y^ the rms- value of velocity fluctuations 
averaged over the convectively unstable layer and time. In 
our calculations Co varies about two orders of magnitude 
from about 0.1 to roughly 15 (see Tables and [2Jl ■ 

Convection efficiency is measured by the Rayleigh 
number 



Ra 



d'gS 



(16) 



where 6 = \7 — Vad is the superadiabaticity, measured 
as the difference between the radiative and the adiabatic 
logarithmic temperature gradients, and Hp the pressure 
scale height, both evaluated in the middle of the unstably 
stratified layer in the non-convecting hydrostatic reference 
solution. The Schwarzschild criterion states that Ra > 
for the convective instability to occur but in reality it must 
exceed some threshold value Rac. Rotation is known to in- 
crease Rac in general fChandrasekhar il961jl and its value 



at different latitudes varies (e.g. Hathaway & Somerville 
[T983). The Rayleigh numbers (2.5 • 10^ and 10*^) used in 
this study are amply supercritical. 

The magnetic field strength is given by the 
Chandrasekhar number 



Ch 



(17) 



where Bq is the magnitude of the imposed field. We apply 
pseudo- vacuum boundary conditions 



Bx — By — — 
dz 



Re 



Rm 



We define the Reynolds numbers as 
utd 

V 

Utd 



(18) 

(19) 
(20) 



The initial stratification is polytropic, described by 
the indices toi, m2, and 7713 for the three layers. We set 
mi — 00, 7712 — 1, and 7713 — 3, respectively, which means 
that the cooling layer is initially isothermal, and the strat- 
ification of the lower overshoot layer resembles that of the 
solar model of Stix l)2UU2|l . For the convective instability 
to occur, the Rayleigh number must be positive, which 
requires that mi < mad — 3/2. The stratification for the 
internal energy e can thus be written as 



e{zo < z < zi) = eo , 
e{zi < z < Z2) = ei + 



'{z-zi) 



(7-l)(m2-Hl) 



(21) 
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e(z2 < Z < Z3) = 62 + 



(7-l)(m3 + l) ■ 
Corresponding equations for the density are 



p{zo < z < zi) 
p{zi < z < Z2) 

P{Z2 < ^ < Z3) 



Po exp 

14 



Pi 

P2 



gjz- zq) 

(7 - l)eo 
9{z- zi) 



1 



(7 - l)(m2 + l)ei 
9{z- Z2) 



(22) 



(7 - l)(m3 + 1)62 



Initially the radiative flux, i^rad — kV6, where k — 
JPX is the thermal conductivity, carries all of the energy 
through the domain. This constraint together with the 
Eqs. ()21|l define the thermal conductivities in each layer 
as 



nij + 1 
rrii + I 



(23) 



where rrii and mj are the polytropic indices of the re- 
spective layers. In the calculations the vertical profile of k 
is kept constant, which with the boundary condition for 
the internal energy, Eq. H12|l . assures that the heat flux 
through the domain is constant at all times. Furthermore, 
K is smoothed at the interfaces of the layers over an inter- 
val of about O.ld by a sine function where the inflection 
point occurs at the gridpoint nearest to the interface. The 
initial state is perturbed with small scale velocity fluctu- 
ations of the order of 10~'^^/dg which are deposited in the 
convectively unstable layer. 

The code used in the calculations is a modified ver- 
sion of that presented in Gaunt & Korpi l|2()()l|l . The 
numerical method is based on sixth order accurate ex- 
plicit spatial discretisation and a third order accurate 
Adams-Bashforth-Moulton predictor-corrector time step- 
ping scheme. The code is parallelised using message pass- 
ing interface (MPI). The calculations were carried out on 
the IBM eServer Cluster 1600 supercomputer hosted by 
CSC Scientific Computing Ltd., in Espoo, Finland. 

3. Summary of the computations, tests, and 
methods 

The parameters used in the calculations are listed in 
Tables n and |2] along with a few diagnostic quantities cal- 
culated from the data. The calculations produced are la- 
belled with the Coriolis number, estimated using the tur- 
bulent velocity from a non-rotating calculation, and lat- 
itude; as the turbulent velocity is changing as function 
of latitude and rotation especially in the rapid rotation 
regime, the actual Coriolis number realised in a calcula- 
tion does not exactly coincide with its label (see the fifth 
column) . The investigated Coriolis numbers vary between 
0.1 and 15; calculations with even higher values would 
have required higher spatial resolution (see below). 

To cut the computational expense some compromise 
concerning the latitudinal coverage had to be made. In 



sets CoOl, Col, Co4, and ColO calculations were made 
at latitudes 0, -7.5°, -15°, -30°, -45°, -60°, -75°, and 
—82.5°, whereas in the sets Co05, Co2, and Co7, only 
at latitudes -7.5°, -30°, -60°, and -82.5°. On the poles, 
the Reynolds stresses should vanish due to symmetry and 
thus the runs there were limited to a few test cases to 
confirm the expected behaviour. 

The numerical resolution was kept fixed at a rather 
moderate 64'^ grid points for the sets CoOl, Co05, Col, and 
Co2. This resolution turned out to be insufficient for the 
rapid rotation cases where strong azimuthal mean fiows 
were generated. Since the spatial resolution in the ver- 
tical direction was already higher than in the horizontal 
directions and the vertical velocity was less affected by 
rotation (compare columns 9 and 10 in Tables ^ and , 
increasing the horizontal resolution moderately was a suf- 
ficient solution. Thus, for the full sets Co7, and ColO, and 
a subset of the Co4 runs a resolution of 96^ x 64 was used. 
Since one calculation in the ColO set still ended in a nu- 
merical instability (ColO-15), we decided not to increase 
rotation further with the expense of even higher compu- 
tational needs in the present study. Comparison with the 
two resolutions can be made within the Co4 set. The dif- 
ferent resolutions are labelled with the suffix 1 (for lower 
resolution) and h (for higher resolution) in Tabled The 
correspondance of the calculated quantities is satisfactory 
between the low and high resolution cases. We feel that 
the results obtained this way are consistent. 

To see how the turbulent dynamics depends on Ra, 
we have performed subsets (latitudes —7.5°, —30°, —60°, 
and -82.5°) of sets Col, Co4, and ColO with a Rayleigh 
number of 10^. The numerical resolution was 128^ for the 
HICol and HICo4 sets and 192^ x 128 for the set HIColO. 
See Table |21 for further details. In what follows, we will 
present results mostly from the Ra = 2.5 • 10^ runs and re- 
port on the differences between the two Rayleigh numbers 
when they are significant. 

3.1. Mean and fluctuating quantities 

When measuring the turbulent properties of a flow one 
must make a clear distinction between the mean and fluc- 
tuating parts of the quantities involved. In this paper we 
use the following definitions: the mean of a quantity qi is 
taken to be the horizontal average 



1 



LxLy 



qidxdy 



and the fluctuation is 



(24) 



(25) 



where qi is the total value of the quantity. Using the 
Reynolds rules, the mean value of a correlation of two 
fluctuating quantities q'^q'j can be written as 



(26) 
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Table 1. Summary of the Ra = 2.5 • 10^ calculations. In all cases investigated Pr = 0.4, 6 = 0.1, Hp = 0.45, and 
eo = 0.3. With these parameters v = 5.96 • 10~^. The full computational domain extends over about 4.2 pressure and 
2.7 density scale heights. The averages are calculated over the unstable layer and time. 
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129 


3.25 


10* 


4.62 
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0.077 


0.042 
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0.050 


30.8 


96 
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132 
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10* 


4.54 


0.079 


0.079 


0.043 


0.042 


0.050 


26.7 


100 
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127 


3.25 


10* 


4.71 


0.076 


0.076 


0.039 
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0.051 


28.8 
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126 


3.25 


10* 


4.75 


0.075 


0.075 


0.039 


0.039 


0.051 


32.7 


95.7 
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128 


3.25 


10* 


4.69 


0.076 


0.076 


0.040 
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0.051 


30.4 


93.9 


C07-7.5 


-7.5° 


173 


9.97 


10* 


5.98 


0.105 


0.103 
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0.060 
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134 


131 
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-30° 


142 


9.97 


10* 


7.26 


0.085 
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0.050 


0.048 
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107 
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-60° 


123 


9.97 


10* 


8.51 


0.074 


0.074 


0.041 


0.039 
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-4.5 


89.9 


Co7-82.5 


-82.5° 


123 


9.97 


10* 


8.60 


0.074 
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0.039 


0.039 


0.048 


1.1 


89.1 


ColO-00 


0° 


321 


2.03 


10® 


4.81 


0.248 


0.191 


0.034 


0.142 


0.112 


-18.3 


238 


ColO-7.5 


-7.5° 


172 


2.03 


10" 


8.65 


0.104 


0.102 


0.055 


0.059 


0.062 


93.9 


131 


ColO-15 


-15° 


148 


2.03 


10^ 


9.94 


0.089 


0.088 


0.047 


0.049 


0.055 


-57.1 


33.8 


ColO-30 


-30° 


136 


2.03 


10*^ 


10.9 


0.081 


0.081 


0.047 


0.044 


0.049 


-46.4 


90.2 


ColO-45 


-45° 


117 


2.03 


10® 


12.7 


0.070 


0.070 


0.041 


0.037 


0.042 


-33.3 


75.6 


ColO-60 


-60° 


112 


2.03 


10® 


13.3 


0.067 


0.067 


0.038 


0.036 


0.041 


-20.0 


76.2 


ColO-75 


-75° 


96 


2.03 


10® 


15.7 


0.058 


0.058 


0.031 


0.029 


0.038 


-13.3 


22.6 


ColO-82.5 


-82.5° 


118 


2.03 


10® 


12.5 


0.071 


0.070 


0.041 


0.038 


0.041 


-25.2 


79.0 



from which one can deduce the rms value of the fluctuating 
quantity as 



Furthermore, the volume average of quantity g, over the 
unstable layer is defined by 



= — ^ — / kli}dz . (28) 
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Table 2. Summary of the Ra = 10^ calculations. In these runs v = 2.98 • 10^"*, other parameters as mdicated in Table 
□ 



Run 


e 


Re 


Ta 


Co 




\ iv 








H[10"*] 


AT/r 


HICol-7.5 


-7.5° 


263 


5.79- 


10" 


0.92 


0.081 


0.078 


0.038 


0.039 


0.056 


5.1 


45.9 


HICol-30 


-30° 


266 


5.79 ■ 


lO'* 


0.92 


0.081 


0.079 


0.040 


0.041 


0.055 


9.8 


46.5 


HlCol-oO 


—60 


269 


5.79 ■ 


10* 


0.91 


0.080 


0.080 


0.041 


0.041 


0.055 


15.6 


46.9 


HlOol-oz.o 


—82.5° 


265 


5.79 ■ 


10* 


0.93 


0.079 


0.079 


0.040 


0.040 


0.054 


11.0 


46.3 


HICo4-7.5 


-7.5° 


301 


1.12 ■ 


lO'' 


3.60 


0.090 


0.090 


0.045 


0.054 


0.055 


-7.8 


52.5 


HICo4-30 


-30° 


267 


1.12 ■ 


10^ 


4.11 


0.080 


0.080 


0.044 


0.046 


0.047 


6.9 


46.6 


HICo4-60 


-60° 


245 


1.12 ■ 


10*^ 


4.55 


0.073 


0.073 


0.039 


0.039 


0.047 


26.7 


42.7 


HICo4-82.5 


-82.5° 


237 


1.12 ■ 


10*^ 


4.71 


0.071 


0.071 


0.037 


0.037 


0.048 


28.8 


41.3 


HIColO-7.5 


-7.5° 


328 


7.03 ■ 


10^ 


8.38 


0.098 


0.098 


0.053 


0.056 


0.059 


94.0 


41.7 


HIColO-30 


-30° 


241 


7.03 ■ 


10^ 


11.4 


0.072 


0.072 


0.042 


0.040 


0.041 


-38.5 


35.1 


HIColO-60 


-60° 


205 


7.03 • 


10'^ 


13.7 


0.061 


0.061 


0.034 


0.032 


0.039 


-13.6 


30.0 


HIColO-82.5 


-82.5° 


201 


7.03 • 


10*^ 


14.0 


0.060 


0.060 


0.032 


0.031 


0.040 


-8.4 


28.4 



Additionally, a time average is defined by 



Convergence of <Q„>, 



1 



(29) 



where AT is the length of the average. The Reynolds rules 
are only approximately valid for time averages. However, 
the longer the averaging interval, the better the approxi- 
mation. In the present study the time averaging intervals 
are long enough for the Reynolds rules to be very nearly 
obeyed. All the calculations were started from the same 
initial convectively unstable state from which it in most 
cases takes about 50 time units for the convection to reach 
a saturation. In order to remove the effects of the initial 
transient we always neglect the first 100 time units from 
the calculation of the time average. For rapid rotation sets 
Co7 and ColO with the lower Rayleigh number, the tran- 
sient is somewhat longer which is taken into account in 
the calculation of the averages. In what follows, the time 
average is applied always in addition to a horizontal or 
volume average. 



0.0001 - 



-0.0001 r 




Fig. 2. Convergence of the Reynolds stress component 
(Qxy) from the run Col-45. The approximate number 
of turnover times in the average of the corresponding 
linestyle is indicated in the plot. 



3.2. Statistics and error estimation 

It is known from previous studies (e.g. Pulkkinen et al. 
119931 Chan I2001|l that the Reynolds stresses are highly 
fluctuating quantities and thus large data sets and as much 
averaging as possible is needed to derive reliable results. 
Chan H2001|l estimates that reasonable convergence of the 
Reynolds stresses is obtained by averaging at least over 
a hundred turnover times. Fig. |21 shows the horizontally- 
averaged Reynolds stress component (Qxy) from the run 
Col-45, averaged over six different time intervals, corre- 
sponding to 4, 19, 38, 60, 81, and 102 convective turnover 
times. The convergence of this quantity is indeed quite 
slow in this case, but the main features of the stress do 
not change significantly if the average is taken over more 
than 40 turnover times. 

We list the integration times (minus the initial tran- 
sient) of the calculations in Tables ^ and |21 from which 



we see that the integration time for the bulk of the runs 
exceed these limits, except for the cases which were termi- 
nated due to numerical problems. In order to restrict the 
computational costs, the Ra = 10^ runs were integrated 
only approximately half the time (in dimensionless units) 
in comparison to the low resolution calculations. 

We estimate the error of the averaged quantities using 
a modified mean error of the mean, defined as 



1 



N 



AT (A^-1)^ 



1/2 



(30) 



where is the number of turnover times in the aver- 
age, giving an estimate of the number of independent re- 
alisations of convection during a calculation, and N the 
number of data points in the average (/). 
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4. The A-effect and turbulent heat transport 

4.1. A-effect 

We begin by introducing the mean field momentum equa- 
tion (the Reynolds equation) in the magnetohydrody- 
namical regime. Firstly, the velocity and magnetic fields 
are separated into the mean and fluctuating parts, U — 
(u) + u' and B = (B) + 6', which are then inserted into 
the momentum equation, Eq. (jS)). Taking the mean and 
applying the Reynolds rules one arrives at 



d(u) 
dt 



-V-{pQ-M) + pg-Vp (31) 
-(J) X (B)-V-(a). 



where the fluctuations of p and p have been ignored. The 
influence of the small-scale turbulence on the large-scale 
velocity comes through the Reynolds stress Qi 



and the Maxwell stress Af^ = ^^^{b'fi'j) tensors, whilst 
the influence of the large-scale magnetic field through the 
Lorentz force (J) x (B), where (J) = P-o^^ ^ (B) is the 
mean electric current. 

Furthermore, by adopting a spherical coordinate sys- 
tem, assuming that the quantities are axisymmetric due 
to which defining the mean values as azimuthal averages is 
meaningful, and neglecting the molecular viscosity which 
is known to be small in stellar convection zones, one can 
derive the equation for the angular momentum balance 



d 

—{ps^n) + \/-{ps^nur, 



s{pu'^u') 
p~\b'^b'))^0, 



(32) 



where s = rsinO is the lever arm, Um = {{ur) , (ue) , 0) 
the meridional circulation, and (u^) — rfl sin 6 describes 
the basic rotation. From this equation it is evident that 
both the meridional circulation and the turbulent stresses 
contribute to the angular momentum transport. 

If the Boussinesq ansatz was adopted for the Reynolds 
stresses, i.e. writing them in the form of a diffusion term 
Qij = ~^t( ^g"'^ + they would act to reduce the 

gradients in the angular velocity rather than as a source 
of differential rotation. This, however, leads to a contra- 
diction between the observed Reynolds stress component 
Qg^ from sunspot proper motions (Ward 119651 more re- 
cently determined by e.g. Pulkkinen & Tuomincn 1998) 
and the one calculated from the observed latitudinal dif- 
ferential rotation. As proposed by Riidiger l)1977|l in the 
form of the theory of the A-effect, the rotational influence 
on a turbulent anisotropic flow may give rise to additional 
stresses proportional to fl itself, capable of driving differ- 
ential rotation rather than diffusing it. Taking this effect 
into account, also the sunspot measurements can be ex- 
plained (e.g. Tuomincn & Riidiger 1 19891 Pulkkinen et al. 

The theory of the A-effect is derived using the mean- 
field approach for hydrodynamics (see e.g. R89 for a 
thorough discussion). In a statistically steady and non- 
magnetic medium, the Reynolds stresses may be expanded 



in increasing powers of spatial derivatives of the mean ve- 
locities 



Qij — ^ijk^k — J^ijkl 



d{uk) 
dxi 



(33) 



where A^fe and Miju are tensors which contain the non- 
diffusive and diffusive parts of the stresses, respectively. 
Higher order derivatives of the mean velocity can be ne- 
glected under the assumption of scale separation (i.e. that 
the mean- field approach is applicable). In the present 
study we confine our interest to the non-diffusive part of 
the Reynolds stresses, but note here that the diffusive part 
has been discussed in length by R89 and more recently in- 
vestigated with analytical turbulence models based on the 
quasi-linear approach (e.g. Kitchatinov et al. I1994b)l . 

An appropriate form of the A-tensor leading to non- 
diffusive expressions can be constructed if the anisotropy 
due to the stratification in the radial direction is taken 
into account (see e.g. R89), resulting in 



/ 

y Ay sin 9 Ah cos ( 
Ay = iyt{V^'>^r,n) + Vir,e,n)) ,and 
Ah = :^tH{r, B, 9) . 




n , where (34) 

(35) 
(36) 



Here is the turbulent viscosity, 9 the colatitude, and the 
coefficients V and H describe the generation of vertical 
(radial) and horizontal (latitudinal) differential rotation. 
In the limit of slow rotation only the fundamental mode 
is expected to be important, whilst for higher rotation 
rates also the H and higher V modes nonlinear in may 
become excited (e.g. KR93). 

The rfl-component of the Reynolds stress tensor does 
not enter to the equation of angular momentum transport 
(Eq. 1321) . but appears in the meridional components of 
the Reynolds equation (Eg. 1^ . thereby being capable of 
influencing the meridional circulation pattern. As noted 
e.g. by R89, such a correlation may arise from the ten- 
sor iliilj, and can therefore be expected to take the form 
Qre sin 6* cos 6*. In analogy with the Eqs. H35|) and H36|l 
we write 



ro — vt AI{r, 9, il) sin 9 cos 9 



(37) 



These coefficients are now obtainable from the calculations 
via the equations 



H = - 



V = - 



M 



{Qxy)v 

vt cos 99 

{Qyz)v 
Vt sin 99 

{Qxz)v 



(38) 
(39) 
(40) 



Vt sin 9 cos 99 

where the minus signs are due to the transformation from 
the coordinate system used in the calculations to the 
spherical ones. 

The mean velocity field can be afftected by the pres- 
ence of a dynamically significant magnetic field by two 
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principal ways, firstly through the large-scale Lorentz 
force iJ,Q^{J) X (B), also known as the Malkus-Proctor ef- 
fect, secondly through the turbulent Maxwell stresses , 
both of which are, in general, thought to reduce the gener- 
ation of differential rotation through the A-effect; thereby 
this effect is referred to as the A-quenching. In some in- 
vestigations (e.g. Jennings & Weiss [TMTjl this quenching 
was assumed to take a similar form to the a-quenching of 
mean- field dynamos, namely ^ cx ^/{l + kB"^). One 
must, however, note that it is usually not the differential 
rotation itself but the A-coefRcients that are quenched in 
this simple algebraic fashion. More detailed investigations 
have been elaborated in the framework of the quasi-linear 
approach (e.g. Kitchatinov et al. 1994b), indicating that 
the behaviour of the turbulent stresses in the magnetic 
regime is complex and cannot be described by such a sim- 
ple quenching formula. In contrast, in the regime of a weak 
magnetic field, the stresses can actually be enhanced, al- 
though in the magnetically dominated regime strong sup- 
pression is predicted (for more details see Sect. l5.7.1|l . This 
effect has been taken into account in some solar mean-field 
dynamo models, and it has been found to play an impor- 
tant role, e.g., for reproducing the torsional oscillations 
(Kiiker et al. lTM^ . minima of magnetic activity reminis- 
cent of the observed solar minima (Kichatinov et al. .l994bk 
Kiiker et al. 11999(1 and longer term variation similar to the 
Gleissberg cycle (Pipin[ 



4.2. Turbulent heat transport 

The mean energy equation in terms of the temperature, 
T, can be obtained from Eq. ^ by applying the Reynolds 
rules similarly as above in the case of the momentum equa- 
tion 



fd{T) 



v(r) 



~V ■ {pcv{u'T'})~pV ■ {u) 
-V(.cV(T» + (r) , (41) 



where (F) contains the contributions from viscous and 
Joule heating and the cooling near the surface. The turbu- 
lent heat transport is described by the correlation tensor 
{u'^T'), which can be expected to be latitude dependent 
due to the influence of rotation on the turbulence. This 
latitude dependence can induce temperature differences 
which are capable of driving meridional flows, which in 
turn may affect the angular momentum transport (as ev- 
ident from Ea. l^ . 

Following the treatise of R89, this tensor can be rep- 
resented by 



(42) 



where Xij is the eddy heat conductivity tensor and Sj = 
— (S/jT — gj/cp) the superadiabatic temperature gradient 
in the direction j . In the present model geometry the hor- 
izontal temperature gradients vanish due to the periodic 
boundaries, and we must thus limit the discussion to the 
meridional components of the tensor. In analogy to the 




Fig. 3. Snapshot of temperature fluctuations and velocity 
vectors from the run HICol-SO. At the top surface of the 
box we show horizontal contours at the z — zi plane and 
at the lower surface at the z = zi plane. The coordinate 
system and the direction of the rotation vector are shown 
in the upper right corner. 



A-coefficients, Eqs. and we write Xij as 



X9r = -XtHV{r,e,^) 



(43) 
(44) 



where Xt is the turbulent heat conductivity for which we 
assume that xt = i^t- 

The correlations between velocity and temperature 
fluctuations derived from the calculations can now be used 
to determine the coefficients VV and HV with the equa- 
tions 



VV 



HV ^ 



V 



PzXt 



PzXt 



(45) 
(46) 



where the sign changes for both quantities due to the 
transformation to spherical coordinates. 

5. Results 

5.1. Convective structures 

To illustrate the general appearance of the structures aris- 
ing in the present calculations we show a typical snapshot 
of temperature fluctuations overplotted with the veloc- 
ity vectors from fully developed convection in a moderate 
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Fig. 4. Horizontal cross-sections of temperature from fully developed convection at the top of the convectively unstable 
layer, z = z\. From left to right, the individual frames are from latitudes 0, —30, —60, and —82.5 degrees. The rotation 
rate increases from top to bottom, corresponding to sets CoOl, Col, and ColO. Here, the x-axis points to the east and 
the y-axis to the north. 




Fig. 5. Vertical cross-section of temperature fluctuations 
and meridional velocity vectors in the xz-plane at y = — 1 
in the run Co7-60. The arrow at the upper right corner 
indicates the direction of the rotation vector. 



rotation case HICol-30 in Fig. |3| As in many previous 
studies we find a convection pattern dominated by broad 
warm upflows and narrow cool downflow plumes. Near 



the top the upflows form cells, separated from each other 
by a network of downflow structures at the cell bound- 
aries. The downflow structures remain coherent over large 
depths, sometimes even extending over the whole convec- 
tively unstable layer. Generally, the horizontal scale of 
the convective pattern is observed to decrease as a func- 
tion of depth. Note that the downflow lanes seen near 
the top tend to form a few well-defined plumes at larger 
depths and the regions of upward flow connect and be- 
come broader. However, the upflows also tend to become 
less well-defined, so that much more small-scale structure 
is seen as opposed to the large and well-defined 'gran- 
ules' near the top. Although this apparently contradicts 
the prediction of the mixing length theory, on average the 
upflows occupy a larger area in the deeper layers in accor- 
dance with the basic mixing length concept. 

The effects of rotation can be studied from Fig.Qlwhere 
we show horizontal contours of temperature from three 
different calculation sets at four latitudes. In the weak ro- 
tation case (set CoOl), shown in the uppermost row, the 
convective cells have more or less angular shapes, which 
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pattern changes only slightly as function of latitude. There 
is a weak tendency of the sizes of the structures to get 
smaller towards the pole, which effect becomes more pro- 
nounced as the rotation is increased (the middle row, set 
Col); the size of the structures clearly decreases as func- 
tion of increasing latitude. The angular appearance of the 
cells also changes due to the strong vortical downflows 
which are generated by the interaction of the converg- 
ing flow with the Coriolis force at the cell edges. In the 
rapid rotation case (the bottom row, set ColO) the convec- 
tive pattern is dominated by rather irregular small-scale 
structures, except for the equatorial case, where the pat- 
tern is totally aligned with the rotation axis, reminiscent 
of the banana cells seen in some global convection mod- 
els (e.g. Brun & Toomre l2002|l . Such an alignment is a 
generic feature of the convective motions in the rapid ro- 
tation regime, as can be seen from Fig. |S1 where a;2-slices 
of temperature fluctuations overplotted with velocity field 
vectors are shown from the run Co7-60. 

5.2. Velocity field characteristics 

In Tables Q] and |2] we have calculated some diagnostics 
of the velocity field. Comparing the total and fluctuat- 
ing rms- velocities (columns 6 and 7, respectively), one 
can note that both tend to decrease towards the poles 
and increase steeply near the equator for rapid rotation. 
A similar trend can be seen in the azimuthal and verti- 
cal velocity field components (columns 9 and 10). This is 
partly due to the fact that the convection is more efficient 
in the equatorial regions compared to the polar ones (see 
Sect. 15. 4() . partly because the increasing rotational influ- 
ence generates azimuthal mean flows, seen as a deviation 
between the total and the fluctuation, which are strongest 
near the equator. 

The degree of horizontal anisotropy, measured by the 
quantity 

Ah = " 2 . 47 

is very small at all latitudes for the slow rotation cases 
(sets CoOl to Col). As rotation is increased. Ah obtains 
positive values near the equator due to the strong az- 
imuthal mean flows generated there. However, the isotropy 
is more or less re-established at southern latitudes 9 > 
30°, also for the rapid rotation cases. 

The vertical component of the velocity is dominating 
over the horizontal ones for the slow rotation cases, which 
means that the preferred direction of the convective mo- 
tions is the radial one, in which stratification and energy 
transport also occur. This anisotropy, measured by the 
quantity 

..^<<>^, (48) 

is observed to decrease as the rotational influence in- 
creases, because the convective structures become more 
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Fig. 6. Volume-averaged kinetic helicity, H, normalised 
by the maximum value, for the sets Col (dashed line) and 
ColO (solid line). For comparison, the dotted line shows 
a — cos 9 profile and the dash-dotted line a — sin 9 profile 
multiplied by 0.65 in order to fit to the data. 

and more aligned with the rotation vector. The sign of 
this quantity is also useful since it should coincide with 
the sign of the radial gradient of the angular velocity, ^ 
(Riidigeril980(). In the present calculations. Ay is almost 
always negative, except for the rapid rotation cases near 
the equator. Thus, the velocity anisotropy generated in 
our calculations reproduces the helioseismic result of pos- 
itive ^ at the equator and at low latitudes and a nega- 
tive one at higher latitudes only for the sets Co4 to ColO. 
However, we expect that this result is slightly affected 
by the fact that the averages were taken only over the 
convectively unstable region. In this case the flow regions 
where the vertical flows turn into horizontally diverging 
flows, mostly occurring in the upper boundary and over- 
shoot layers, were excluded from the calculation, resulting 
in the overestimation of Uz over Uy . 

5.3. Kinetic helicity 

Column 11 in Tables Q and |21 shows the volume averaged 
kinetic helicity, Ti., defined by 

n = {(y X u') ■ u') . (49) 

This quantity is interesting in that it is believed to play 
a major role in the dynamo theory, providing the so- 
called a-effect, which describes the generation of large- 
scale poloidal magnetic fields due to the inductive action 
of helical turbulence. In the simplest approximation, mak- 
ing use of the first order smoothing and assumption of 
isotropy, a ~ — ir'Ti, where r' is the correlation time as- 
sociated with the turbulence. It has been noted that the 
a computed from magnetoconvection calculations qualita- 
tively agrees with the above approximation (Ossendrijver 
et al. 120011 12002|l . Furthermore, it is possible to determine 
the sign and latitudinal distribution of the magnitude of 
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Fig. 7. Mean energy fluxes as functions of depth from the 
run Col-30, normalised by the input flux, see Sect. l5.4l for 
definitions. 



Ti. in the surface layers of the Sun from helioseismology 
(Duvall & Gizon i2000. see also Riidiger et al. 1999 who 
studied numerically the determination of H from surface 
flows). The aforementioned studies indicate that Ti is max- 
imal at the poles with a nearly cos 0-latitude dependence. 
The sign of H is noted to be positive in the bulk of the 
convection zone in the southern hemisphere and negative 
in the northern hemisphere. 

We find that on the low rotation calculations H shows 
several sign changes as function of depth so that the vol- 
ume average is small and does not show a clear latitudinal 
trend. When rotation is moderate (sets Col to Co4) there 
is a clear tendency for the helicity to be positive in the 
convectively unstable region and to increase towards the 
pole (bar the Co4 set with Ra 2.5-10^), in accordance with 
the aforementioned studies. As function of depth, the posi- 
tive values peak near the top of the convection zone while 



a region of smaller negative helicity resides in the over- 
shoot layer. This region of negative Ti, grows as rotation 
is increased so that for the strong rotation cases, only a 
narrow region of positive helicity remains in the surface 
layers. Near the equator, this positive peak dominates in 
the volume average of Ti, but for latitudes 8 < —15° the 
sign is negative which behaviour carriers all the way to the 
pole, however, with monotonically decreasing absolute val- 
ues (see Figure E}. Taken that the close relation between 
the helicity and the a-effect persists in the rapid rota- 
tion regime, this would indicate that for strong enough 
rotation, the sign of the a-effect changes from negative 
to positive and simultaneously decreases in magnitude at 
higher latitudes, whilst larger negative values occur in the 
equatorial regime. Interpreted in terms of the classical dy- 
namo theory, in the regime of slow and moderate rotation, 
the poloidal field generation would be more efficient near 
the poles with equatorward migration of activity phenom- 
ena, when combined with the observed internal rotation 
law. In the rapid rotation regime, in contrast, the a-effect 
would be more efficient in the equatorial regions, and the 
migration would be equatorward at low latitude regions 
and poleward at higher latitudes. However, the propor- 
tionality between Ti. and the a-effect is only approximate, 
and can yield incorrect results for realistic circumstances. 
Thus, calculations with magnetic fields are needed in or- 
der to determine whether the behaviour suggested by the 
helicity really occurs in the rapid rotation regime. 

5.4. Convection efficiency 

To show the energy transport in a typical calculation, we 
plot the depth dependence of the horizontally averaged 
energy fluxes in Fig. for the run Col-30. We define 

i^rad = , (50) 

oz 

Fconv = -j{pe'u'^) , (51) 
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F ■ 

^ vise 



(52) 
(53) 



as the radiative, convective, kinetic, and viscous fluxes. 
The minus signs arise because we define positive fiux to 
be directed outwards. Most of the energy is transported 
by radiative diffusion, whilst the convection transports 
roughly one fifth of the total energy flux through the con- 
vection zone. The fiux of kinetic energy is directed down- 
wards, and amounts to a few per cent of the total flux. The 
viscous flux is in all cases negligible. The change of sign 
of i^conv in the overshoot layer is due to the fact the nega- 
tive temperature fluctuations entering the stable layer flnd 
themselves as positive fluctuations in the new background 
stratification. We find that the convective flux decreases 
as rotation is increased so that i^conv in the ColO set is 
roughly half of that in the CoOl set. 

The rotational influence on the convection efficiency 
can be estimated quantitatively by the Nusselt number, 
(see e.g. Brandenburg et al. il99n|l 



Nu = 



Fad 



F. 



(0) 



ad 



F, 



ad 



where Fw 



F, 



ad 



F,. 



F, 



kin 



(54) 



is the total fiux 



through the convectively unstable layer, F^d — '*(§f) 



ad 



the radiative flux if the layer is adiabatic, and F 



rad 



Az 



the radiative fiux if the temperature stratification in the 
unstable layer is linear. The Nusselt number gives the ratio 
of the total heat flux and the heat fiux which would be 
transported by conduction alone. Thus Nu must exceed 
unity for convection to be present. 

The resulting Nusselt numbers as functions of rotation, 
measured by Ta and Co, at the latitudes investigated are 
shown in Fig. |H1 The slow rotation case CoOl does not 
show any definite latitudinal variation, which is indeed 
the expected behaviour. The variation seen within the set 
CoOl can then be used as an error estimate for the other 
sets. When the rotation is stronger, the efficiency of con- 
vection decreases as a function of latitude in a more or 
less monotonic fashion. Convection at equatorial regions 
seems to be most vigorous and for 9 < —45° the efficiency 
does not change appreciably as the latitude increases. The 
trend that convection at polar regions would be more 
vigorous than at some intermediate latitudes noted by 
Gilman l)1977|l and Pulkkinen et al. (1993) is not observed 
here. Furthermore, we do not find a regime where Nu in- 
creases with increasing Ta (Hathaway & Somerville'1 983(l . 
For the Ra =10^ calculations the qualitative results are 
very similar, the distinction being the roughly 50% larger 
values of Nu. 

5.5. Overshooting depth 

Convective overshooting is thought to be important for ex- 
ample in the context of stellar magnetism providing means 
to store magnetic flux below the convection zone by resist- 
ing the buoyancy effect. In our calculations, indeed, the 



downflows very often overshoot into the stable layer be- 
low, where they decelerate due to the buoyancy force. We 
deflne the depth of overshooting as 



Z2 — Zos, 



(55) 



where Zqs is the depth where the horizontally averaged ki- 
netic energy flux falls below a threshold value of e^^ of 
its value at the bottom of the convectively unstable layer. 
In order to facilitate a comparison to a real star we nor- 
malise the overshooting depth by the pressure scale height, 
Hp2 — 62(7 — l)/3, at the bottom of the convectively un- 
stable layer. 

The calculated depths for the Ra 2.5 • 10^ runs are 
shown in the upper panel of Fig. El as functions of lati- 
tude and Coriolis number. As a general feature we flnd 
that the overshooting decreases as a function of increas- 
ing rotation at a given latitude. Such behaviour has ear- 
lier been reported by Brummell et al. (2002) and Ziegler 
& Riidiger l)2()()3|l . For the slowest rotation investigated 
(set CoOl) the overshooting extends rather deep, down to 
0.4_ffp2, and the latitude dependence is very weak as ex- 
pected. For moderate rotation (sets Co05 to Co2) there is 
a clear trend of decreasing overshoot depth with latitude. 
However, as rotation is increased further, there is an op- 
posite trend. At the poles, dos seems to stay more or less 
constant, whilst it continues to decrease at the equator. 
This may be interpreted to be due to the stronger influ- 
ence of the Coriolis force on the vertical downflows at the 
equatorial regions, where they tend to align with the rota- 
tion vector. The lower panel of Fig.lHlshows dos as function 
of the Chandrasekhar number for the magnetic runs (see 
Sect. 15. 7'!^ for a detailed description of the runs). We find 
that the overshooting depth is reduced monotonically with 
increasing imposed fleld strength for the present range of 
parameters where the ratio of the gas to magnetic pressure 
in the bottom of the convection zone spans from « 2.6-10^ 
to 260. Increase of dos for small magnetic field strengths, 
noted by Ziegler & Riidiger (Ii2003j) . is not observed. 

For the most rapid rotation, dos varies between about 
0.177p2 at the equator to about 0.2_ffp2 at high latitudes. 
This value is still larger than the estimates from helio- 
seismology (e.g. Monteiro et al. 1994), but the rotational 
influence in our study is most probably still weaker than in 
the deeper layers of the solar convection zone. One must, 
however, bear in mind that the normalised energy flux, 
Ftot/pCg, where Cs is the sound speed, at the bottom of 
the convectively unstable layer is roughly eight orders of 
magnitude larger than in the case of the sun. Furthermore, 

for non-rotating convection, the penetration depth scales 

3/2 

with Uz and the overshooting depth with the square root 
of the input flux (e.g. Zahn ,1991,. Hurlburt et al. HM^ 
Singh et al. I1998|l . The difference between penetration 
and overshooting is that the former extends the superadi- 
abatic region and the latter does not. In the present cal- 
culations the thermal stratiflcation does not appreciably 
change due to the overshooting and thus the latter scaling 
relation should be applicable. However, in our results no 



14 



Kapyla et al.: Reynolds stresses and turbulent heat transport 



0.40 



0.35 - 



0.30 -\ 



^ 0.25 - 



0.20 



0.1 5 - 




90 105 



0.22 - 



0.20 - 



I 0.1 8 - 



0.1 6 - 



0.1 4 - 




Fig. 9. Upper panel: overshooting depth, dos, as function 
of latitude for the calculations with Ra = 2.5-10^. Different 
linestyles are used for each set of calculations, differing by 
Coriolis number, as indicated in the figure. Lower panel: 
rfos, as function of the Chandrasekhar number. The up- 
per a:-axis denotes the corresponding runs presented in 
Table H 



such scaling exists, since the input flux in the Ra=10^ cal- 
culations is half the one used in the Ra=2.5 • 10^ case but 
the overshooting depth decreases only by a factor of ten 
to fifteen per cent for the Col and Co4 sets and actually 
increases for the ColO set (not shown). A possible expla- 
nation of this behaviour is that the mixing length theory, 
assumptions of which were used in the derivation of the 
scaling relation, breaks down when rotation is present. On 
the other hand, it has to be noted that our description of 
the overshoot region is not very realistic, the thermal dif- 
fusivity profile being essentially a smoothed step function, 
which quite probably acts to hinder the efficiency of over- 
shooting in comparison to a smoothly varying profile. 



5.6. Mean flows 

The generation of differential rotation by the presence of 
anisotropic turbulence and uniform rotation is accounted 
for in the theory of the A-effect (e.g. R89) by the interac- 
tion of the anisotropic eddies with the Coriolis force. In the 
rectangular domain used in this study, this effect works to 
generate mean flows, of which of particular interest are the 
latitudinal and azimuthal ones which are indicative of the 
differential rotation in a full spherical domain. However, 
one must bear in mind that the boundary conditions in a 
spherical domain would permit a variety of other effects 
which affect the mean flows, such as horizontal tempera- 
ture gradients, which in the present geometry could only 
be imposed and cannot arise naturally. 

To investigate the depth dependence of the mean flows 
we average the horizontal velocities over a time period 
designated in Tableland over the horizontal directions. 
This procedure yields {ux){z) and {uy){z) for the latitudi- 
nal and azimuthal mean velocites, respectively, which are 
functions of z alone. The vertical mean velocity {uz){z), 
calculated in this way, is very small in comparison to the 
horizontal ones due to the boundary conditions which do 
not permit a net vertical momentum flux. The fact that 
{uz) still has a small finite value is due to the different 
filling factors of up- and downflows. Fig. 1101 shows typical 
examples of these flows plotted at different latitudes for 
slow (CoOl), moderate (Col), and rapid rotation (ColO). 

The predominant feature of the latitudinal flow for 
slow and intermediate rotation is the mostly negative 
(poleward^) flow in the bulk of the convection zone and the 
positive (equatorward) flow in the overshoot layer. This 
trend is strongest for the set Co2, and when rotation in- 
creases, the latitudinal flow starts to diminish and the di- 
rection reverses for the most rapid rotation. For the rapid 
rotation cases (Co7 and ColO) there is a strong equator- 
ward flow confined to a narrow layer near the surface and 
a much weaker poleward flow in the bulk of the convection 
zone. In the overshoot layer only a very weak equatorward 
flow is present. 

The mean latitudinal flows in the slow rotation regime 
agree with the results of Chan (2001) and Brummell et al. 
lfl998,,) . When rotation is stronger, the qualitative charac- 
ter of Chan's results stays the same, in disagreement with 
ours, whereas in the study of Brummell et al. fl998| a 
sign change near the top of the domain is reported (their 
Fig. 6). The increasing magnitude of the latitudinal flow 
towards the equator noted by Brummell et al. (T998| is 
also present in our results for moderate rotation. However, 
the magnitude of the latitudinal flow in our study seems 
to grow for increasing rotation, a feature not present in 
the two aforementioned studies. 

The azimuthal mean flow shows generally very simi- 
lar behaviour, with increasingly negative values (westward 
motion) near the surface and smaller positive values (east- 
ward) in the overshoot layer for slow and moderate rota- 



Note that the box lies on the southern hemisphere 




Fig. 10. Latitudinal (left) and azimuthal (right) mean flows, {ux) and (wj,), respectively, plotted against depth from 
sets CoOl (top, slow rotation). Col (middle, moderate rotation), and ColO (bottom, rapid rotation). The different 
linestyles correspond to different latitudes, which are denoted in the middle figure. On the lowest graph on the right, 
only a part of the equatorial azimuthal velocity is shown due to its large magnitude. 



tion (up to the set Co2). The magnitude of the azimuthal 
flow decreases more or less linearly as function of depth, 
and also as function of increasing latitude. The latter ef- 
fect is reminiscent of the Taylor- Proudman theorem (e.g. 
Chandrasekhar I1961|l which states that motions cannot 
vary in the direction of the rotation axis leading to cylin- 
der shaped contours of the angular velocity. The former 
effect was discussed by Chan (2Q01 ). who showed that the 
dominant term in the equation for the mean azimuthal 
flow at the equator is the horizontal component of the 
Coriolis force, Jq — 2ilcos8w2, which produces a shear 
of approximate magnitude of 



2VL , 



(56) 



where the difference A is taken over the unstable layer. 
Our results follow this relation for the slow and moderate 
rotation cases (sets CoOl to Co2), but this relation fails for 
the rapid rotation regime. However, a more general way of 
explaining the mean flows in the present calculations can 
be easily obtained by considering the Reynolds stresses, 
see Sect. IS.G.TI below. 

Just as the latitudinal flow, also the azimuthal mean 
flow reverses direction in the bulk of the convection zone 
for the most rapid rotation (Co7 and ColO). There is also 



a moderate westward flow in the bulk of the convection 
zone and a layer of negative shear near the surface. A 
narrow layer of retrograde motion exists in the overshoot 
region. This behaviour is very similar to the investigations 
by Brummell et al. l(TM5jl and Chan 



5.6.1. Generation of the mean flows by the Reynolds 
stresses 

Next we study the horizontal momentum balance in the 
model by deriving the equations for the latitudinal, az- 
imuthal, and vertical momentum balance, which are ob- 
tained by taking a horizontal average of the momentum 
equation, Eq. J^Jl 



d_ 

dz 
d_ 

dz 



{{pUxU;,)) -2n(pUy)smQ ^ , (57) 
{{puyu^)) + 2n{pux) sine ^ , (58) 



-^^{{puD) - 2 n{puy) cosQ - ^ + {p)g - 0. (59) 

In the azimuthal equation we have dropped the Coriolis 
term involving the mean vertical momentum (puz), which 
vanishes due to the boundary conditions. Furthermore, 
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Fig. 11. Reynolds stress components {Qxy) /vt, {Qyz) /vt, and {Qxz) /vt as functions of depth for the sets CoOl (top), 
Col (middle), and ColO (bottom), respectively. Different latitudes are represented by the linestyles denoted in the 
middle figure. 



applying the Reynolds rules, the remaining terms of the 
form (uiUj) can be expanded using Eq. H26|l . In this ex- 
pansion terms with (uz) appear, but these are very small 
in comparison to the other mean velocities and can be ne- 
glected. The remaining terms of Eqs. (|57|l to H59() can then 
be used to solve for the mean horizontal velocities 



1 



d 



{{P){«)), 



2n{p) sinQdz 
1 d 

i{p){K'^'z)) > 



2n{p) sineaz 
1 



(60) 
(61) 



2n{p) cose 



oz oz 



At the equator, the latitudinal mean flow vanishes due 
to symmetry, and Eq. (|60|) has a singularity. However, 
Eq. H62|) . which includes contribution of the diagonal 
Reynolds stress component {u'^) in the pressure field, can 
be used instead. We find that the vertical derivatives of 
the Reynolds stress components {u'yu'^) and {u'^u'^) explain 
the generated horizontal mean flows in all non-equatorial 
models. At the equator, Eq. (|62|l reproduces the observed 
zonal mean flow, and is applicable at all latitudes investi- 



gated. Thus, the observed mean flows in the models arise 
due to the non-vanishing transport of momentum driven 
by the Reynolds stresses as expected. 



5. 7. Reynolds stresses and the A-effect 

In this section we discuss the spatial distribution of the 
Reynolds stress components generated in the calculations, 
and their dependence on rotation and large-scale magnetic 
fields. To draw conclusions of the effects of the generated 
stresses on the angular momentum transport and genera- 
tion of differential rotation in a global spherical geometry, 
we parametrise the stresses in terms of the A-effect coef- 
ficients. This is also interesting because the calculations 
presented here are not restricted by the assumptions of 
FOSA, due to which the results may be used to test the 
validity of this approximation. It has to be born in mind, 
however, that our local calculations also have restrictions, 
including the effects of periodic horizontal boundaries, 
limited number of pressure and density scale heights in- 
cluded, and the small Reynolds and Rayleigh numbers, 
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due to which the turbulence modelled is relatively weak 
and the convection rather inefficient. 

5.7.1. The stresses due to the generated mean flows 

As discussed in the previous section, mean flows exhibiting 
radial gradients naturally arise in the calculations due to 
the non-zero Reynolds stresses (they also vary as a func- 
tion of latitude measured from different calculation sets, 
but due to the periodic horizontal boundaries the latitu- 
dinal gradients vanish in the calculations). In this case, 
in the series expansion of the Reynolds stresses, Eq. 1)33(1 . 
the diffusive contribution will become non-zero as well, 
due to which it should be separated before calculating 
the A-coefficients, for which we have adopted the follow- 
ing procedure. Using the radial profiles of the latitudinal 
and azimuthal mean flows generated in the calculation 
of interest as an additional imposed velocity field com- 
ponent stationary in time (i.e. writing the velocity field 
as {ux + u^^^^" (z) ,Uy + u^°^'^{z),Uz) and adding advec- 
tion terms of types uf^'^'^'^duj/dxi and u^duf^'^^^ / dz to the 
equations) , we repeat a few calculations switching off rota- 
tion, and calculate the generated Reynolds stresses. These 
stresses now correspond to the diffusive part in the se- 
ries expansion, denoted here by Q^^''''^' . Subtracting these 
stresses from the ones generated in the full calculations, 
we obtain the non-diffusive contribution 

of = Qg°*) _ Q(f-^) . (63) 

We chose to apply this procedure only to a subset of the 
full calculations, namely for the sets Col and ColO, cor- 
responding to moderate and rapid rotation, respectively, 
and at the latitudes of 0°, -7.5°, -30°, and -60°; the 
results are summarised in Table |3| 

We note that the stresses generated by the shear flows 
are largest in the cases where the imposed flows were taken 
from the equatorial or near-equatorial cases. The stress 
component Q^xy^^^'' is small in comparison to Q^xy'^^ for 
all latitudes in the Col set (note that the stress compo- 
nents Qxy and Qxz should vanish at the equator due to 
symmetry). In the rapid rotation case Qxy ^ is approx- 
imately one flfth of the total stress near the equator and 
about one eight for latitude —30°. For the other two lat- 
itudes the stresses are statistically consistent with zero. 
Since the contribution of the shear is of the same sign as 
the total, we note that it acts as to diminish the contribu- 
tion of the A-effect. We thus conclude that the horizontal 
A-effect derived from the total stress is overestimated by 
maximally 20 per cent near the equator for the most rapid 
rotation cases. 

The stress component Qy^ shows an opposite trend: 
^(s^hoar) general of opposite sign to Qyz"*'' which indi- 
cates an underestimation of the vertical A-effect. This un- 
derestimation is largest near the equator due to the large 
azimuthal mean flows there (see Figure [T^ . Noteworthy 
is that the shear contribution is almost equal to the A 
part at the equator for the ColO case but for the latitude 



Table 3. Comparison of the volume averaged Reynolds 
stresses generated by rotation and dynamically generated 
mean flows versus stresses generated by mean flows alone 
(in units of [10-^]). 
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-9.39 
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0.04 


ColO-30 


8.68 
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ColO-30-shear 


1.00 


0.05 


-0.58 


ColO-60 


2.62 


-0.28 


4.02 


ColO-60-shear 


1.74 


0.30 


-0.15 



—7.5° the shear contribution is already insignificant. In 
conclusion, the vertical A-effect estimated from (5^2°*^ is 
in general underestimated near the equator possibly by a 
factor of two to three but that this effect is significantly 
weaker for all higher latitudes, being maximally a few tens 
of per cents in the moderate rotation cases. The stress 
component qS'"^''' is insignificantly small in comparison 
to (5^2°*^ for all non-equatorial cases studied. 

5.7.2. Radial distribution 

We continue by discussing the vertical (radial) profiles 
of the Reynolds stresses as the Taylor number is varied. 
These are shown for the horizontally averaged off-diagonal 
Reynolds stress components, normalised by an estimate of 
the turbulent viscosity vt = ^Ufd, in Fig. for slow (set 
CoOl), moderate (Col), and rapid rotation (ColO). In or- 
der to visualise the latitudinal variation of the stresses, 
the data from the two latter sets are shown in a polar plot 
in Fig. nn 

The horizontal component, (Qxy), shows several sign 
changes as function of depth and no clear trends are ap- 
parent for slow rotation (sets CoOl to Co05). For moder- 
ate rotation (Col and Co2) larger negative values tend to 
occur in the convection zone near the equator. For more 
rapid rotation, this latitudinal trend is more pronounced, 
and as a function of depth, a positive peak near the up- 
per boundary develops. The magnitude of this feature in- 
creases with increasing rotation. A smaller positive peak 
develops in the overshoot region for the rapid rotation 
cases. The uppermost panels of Fig.ll2lshow more clearly 
the tendency of (Qxy) to obtain its largest value near the 
equator. The maxima of (Qxy) also tends to occur in the 
regions where the vertical flow is deflected to the horizon- 
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Fig. 13. The volume averaged Reynolds stress compo- 
nents (left), {Qrrt,) (middle), and (Qro) (right) in 
units of [10^^], from the Ra = 2.5 • 10^ calculations, 
normalised by vt^. The errorbars are calculated from 
Eq. I|3U|) . For comparison, the dashed lines show the re- 
sults of KR93 for the corresponding Coriolis number. 

tal directions, namely in the overshoot layer immediately 
below the convection zone and near the upper boundary. 

The vertical Reynolds stress component, {Qyz), shows 
a distinct trend of being mostly positive in the bulk of 
the convection zone for the slow and moderate rotation. 
Considerably weaker positive values of the stress occur in 
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the overshoot region for these cases. The magnitude of 
this component increases up to the set Co2, after which 
the trend begins to change and for the set Co4 this quan- 
tity is more or less zero in the bulk of the convection zone 
near the equator (the positive peak in the overshoot layer 
remains) and the overall magnitude at other latitudes de- 
creases. Similar trend continues in the set Co7. There is 
a sign change in the bulk (but not in the stable region) 
for the most rapid rotation (set ColO). As the latitude 
increases, the negative dip near the top of the domain be- 
comes narrower. 

The meridional component of the stress, {Qxz), is neg- 
ative throughout the domain for the slow rotation cases 
CoOl to Col. The magnitude of this quantity is less in the 
overshoot region than in the bulk of the convection zone. 
For the set Co2 the same trend continues on the lower lat- 
itudes but there is a sign change at higher latitudes. With 
even more rapid rotation, (Qxz) still exhibits a negative 
dip in the unstable region, but in addition a comparable 
positive peak appears in the overshoot region. 

In comparison, the depth dependencies of the Reynolds 
stresses are in many cases at variance with those obtained 
by Chan {2001). The accordance seems to be the best for 
the slow rotation cases (Co < 1) and for the components 
{Qyz) and (Qxz)- The sign changes observed when the 
rotation becomes large enough are seen to occur in the re- 
sults of Chan H20U1|I only for (Qxy) , whereas in our case all 
components exhibit this behaviour. The differences here 
also reflect the variance seen in the mean velocities. 

5.7,3, Latitudinal distribution 

Figs.ElandEldepict the volume averages of the Reynolds 
stresses normalised by I't^ as functions of latitude for all 



Kapyla et al.: Reynolds stresses and turbulent heat transport 



19 



Co Co 
0.5 - ? 4 7 10 0,5 1 2 4 7 10 




4.0 4.5 5.0 5.5 6.0 4.0 4.5 5.0 5.5 6.0 

og,QTa log,QTa 



Fig. 15. CoefScients H (left) and V (right) from different 
latitudes as function of the Taylor number. 

sets of calculations. We have adopted spherical coordi- 
nates in order to facilitate comparisons with other stud- 
ies. Here we give an error estimate using a modified mean 
error of the mean, Eq. (|30|) . For comparison, we show the 
results of KR93 for the corresponding Coriolis number for 
each set. 

We begin our analysis from the horizontal component 
{Qe<p)v (shown in the left columns of Figs. [T^ and . 
which shows a lot of scatter and the errorbars are large 
for the slow rotation cases Co < 1; it seems that no real ef- 
fect is present for these cases. As the rotation is increased, 
a trend starts to appear: the sign of the stress is negative 
at all latitudes, and the maximal magnitude is obtained at 

o 

the latitude of = —7.5 accompanied by a rapid increase 
near the equator and a slower decrease when moving to- 
wards the poles; a series expansion in increasing powers 
of sin^ 6 traditionally used in describing the latitudinal 
distribution (e.g. R89) fits poorly to these profiles (not 
shown) . This behaviour is very similar to the one reported 
by Pulkkinen et al. H1993|l for Co « 5 in our notation, and 
visible also in the results of Chan H2001|l . 

For this stress component observations from the Sun 
exist (e.g. Ward lTMHl VirtanenCnHni). On the sun {Qg(),)v 
is negative on the southern hemisphere with a more or 
less a linear variation as a function of latitude (Pulkkinen 
& Tuominen 11998(1 . The negative sign indicates that 
there is an angular momentum flux towards the equator. 
Qualitative correspondence with the observations is fairly 
good when Co > 2. 

A comparison with KR93 shows that, for the moderate 
rotation cases, the sign of {Qg(j,)v is the same in both stud- 
ies, but that the latitudinal distribution is much more con- 
centrated towards the equator in the present results. Also 
the magnitude of the horizontal A-effect, H = ^'3g.») v 

^ ' vt COS 9\l ' 

seems to be larger in our case. We find that the abso- 
lute value of {Qe<i})v increases up to the set Co7, but the 
growth seems to have saturated for the set ColO, where the 
magnitude is comparable to the Co7 set. What this implies 
is that H itself is decreasing as function of Q. for the most 
rapid rotation cases. This is an indication of rotational 
quenching of the horizontal A-effect, see the left panel of 
Fig. El which is noticiable when Co > 4. Recent obser- 
vational studies of rapidly rotating stars seem to support 



the view that generation of horizontal differential rotation 
is strongly suppressed when rotation is large enough (e.g. 
Reiners & Schmitt l20n3b|l . 

The stress component {Qr4>)v (the middle columns of 
Figs. El and I14|l describes the radial flux of angular mo- 
mentum, and is thus vital in generating radial differential 
rotation. We find that {Qr(f>)v is mostly negative, corre- 
sponding to inward transport of angular momentum, for 
all sets except the case CoOl, where no clear trend is seen 
and the two most rapid rotation cases where a change of 
sign occurs. The magnitudes at the midlatitudes tend to 
be the largest for all the sets up to Co4. As the rotation 
is increased further, the sign changes, starting from the 
equatorial regions so that for the set ColO, {Qr4,)v is pos- 
itive at all latitudes, and obtains the maximum near the 
equator. The qualitative picture of {Qr4>)v for the sets Col 
to Co4 is in accord with Pulkkinen et al. H1993|l and Chan 
11200 1(1 . whereas on the rapid rotation end we seem to have 
come to a regime not yet encountered in the other stud- 
ies. A comparison to the results of KR93 seems to support 
this view as well. The relatively strong and positive verti- 
cal A-effect of KR93 for slow rotation is not present in our 
results. When rotation is intermediate (sets Co2 and Co4), 
the correspondence is rather good, but the sign change at 
high rotation is not seen in the study of KR93. We note 
here that absolute value of {Qr<t,)v does not change ap- 
preciably from the set Co2 to ColO, which means that the 
vertical A-effect V — ^'^'■'^In is also subject to strong ro- 

ut sin ^ 

tational quenching, and that this effect is much more clear 
for V than for H (see the right panel of Fig. I15|l . 

Finally, the Reynolds stress capable of driving merid- 
ional ffow, {Qre)v, does not show as clear trends as the 
two other components. For slow rotation (sets CoOl to 
Col) this component is negative at all latitudes. For mod- 
erate rotation (set Co2), the values near the pole start to 
get positive. Towards the rapid rotation end, there seems 
to be a negative dip at low latitudes (8 ~ —15° ... — 30°) 
and a positive one at high latitudes (6 ~ —60°). This 
qualitative trend is again similar to what Pulkkinen et al. 
(ITO^ find in their study. 

Our main conclusions for the Reynolds stress compo- 
nents relevant for the generation of differential rotation 
indicate that the angular momentum is transported to- 
wards the equator when rotation is significant (cases Col 
and up) and radially inwards for slow and moderate rota- 
tion (sets Co05 to Co4), and outwards for rapid rotation 
(sets Co7 and ColO). The efficiency of the angular mo- 
mentum transport is largest at mid-latitudes for slow and 
moderate rotation and near the equator for the rapid ro- 
tation cases. 

The observed rotational quenching of the A-coefficients 
is of great interest in the context of the dynamo theory, 
where the generation of the toroidal magnetic field occurs 
via the f2-effect, dependent on the absolute differential ro- 
tation Afl. According to our results the generators of the 
fl-eSect may be significantly reduced, the quenching of the 
V coefficients being stronger and taking place with lower 
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Fig. 16. The volume averages of the rms-velocities (left) 
and rms-magnetic fields (right) in different directions. 

rotation rates than for the H coefficients. This implies that 
the horizontal differential rotation would dominate over 
the radial one in stars with moderate rotation rates, but 
that even the horizontal differential rotation would dimin- 
ish significantly in the rapid rotation regime. Therefore, 
the dynamos operating in rapidly rotating stars can be 
expected to be of a^-type. One must, however, note that 
the a-effect is probably, at least at the poles, subject to 
rotational quenching as well (Ossendrijver et al. l2UUf|l . It 
remains to be seen how the reverse of the sign and latitude 
distribution of the helicity noted in Sect. 15.31 in the rapid 
rotation regime changes this picture. 

5.7.4. The effect of an imposed azimuthal magnetic 
field 

We investigated the behaviour of the total turbulent 
stresses, defined as 



T,, 



(64) 



where Qij and Mij 



are the corresponding Reynolds and 
Maxwell stress tensor components, in the magnetic regime 
by a set of calculations (listed in Table 0)) in which we im- 
pose a large-scale azimuthal magnetic field B = (0, Bq — 
const, 0) of varying strength on top of a convection pattern 
at a saturated state. This choice seems plausible since the 
dominating component of the solar mean magnetic field 
is the toroidal one. Since the imposed field has no gradi- 
ents, the large-scale Lorentz-force vanishes in the calcula- 
tions, whilst the Maxwell stresses are non-zero due to the 
turbulent magnetic field component generated when the 
convective turbulence is acting on the mean field. We de- 
fine the modified A-coefRcients including the effect of the 
Maxwell stresses as 



Htot — 


{Txy)v 

vt cos 9Vl ' 


(65) 


Vtot = 


{Tyz)v 

vt sin Q£t ' 


(66) 


Mtot = 


{Txz)v 

Vt sin 9 cos OVl ' 


(67) 



and list the values obtained from the calculations in the 
columns 4, 5, and 6 of Tabled 
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Fig. 17. Vertical profiles of the horizontally averaged to- 
tal turbulent stresses (Ty) normalised by the turbulent 
viscosity from the magnetic runs MO to M4 transformed 
into spherical coordinates. 



As a base calculation for this study we took the run 
Co2-60, where all the Reynolds stress components were 
appreciably large. The maximum Chandrasekhar number 
was chosen to be 10^ (run M4) which corresponds to a 
situation where the magnetic energy is roughly four times 
the kinetic energy of turbulence (Tabled) column 8) in the 
saturated state. This calculation can be thought to repre- 
sent the superequipartition case, whereas M3 corresponds 
to equipartition. 

When a magnetic field is imposed on top of the con- 
vection pattern, the flow undergoes a transition stage of 
a few tens of time units, after which it relaxes to another 
stable state; in measuring the quantities listed in Table 0] 
we have neglected this transition phase. Due to the con- 
vective motions acting on the magnetic field, a fluctuating 
magnetic field is generated. The strength of the fluctuat- 
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Run 


Ch 


Rm 




Kct[10-^] 


Mtot[10-2] 


Ekin [10-3] 


-Emag / -Ekin 


Ut 


bt 


Ah [10-3] 


Av [10-3] 


AT/r 


MO 







4.38 


-3.67 


-5.80 


2.03 




0.085 




1.6 


-12.1 


108 


Ml 


1 


141 


-0.01 


-3.57 


-2.41 
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0.22 


0.084 


0.009 
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57.8 
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Fig. 18. A-coefHcients i?tot, Vtot, and Mtot as functions 
of the normalised field strength, (3m = ^/JIopu^- The 
quenching functions derived by Kitchatinov et al. ( 1994bfl 
are overplotted with the dotted linesymbols (see Sect. 
15. 7.41 for explanations). 



ing field (measured by the rms-alue listed in column 10) 
grows roughly exponentially up to M3, after which the 
growth seems to become quenched. For the runs Ml to 
M3 the fluctuating field is dominating over the large-scale 
field, although the fraction of small-scale energy over the 
large-scale one is decreasing nearly exponentially when the 
strength of the imposed field is increased. For the run M4 
the large-scale field is dominating over the fluctuating one 
and also tends to elongate convective structures along its 
direction. The total and fluctuating velocity fields, mea- 
sured by the total kinetic energy and the rms-value listed 
in columns 7 and 9, respectively, become suppressed by 
the magnetic field, as expected; from Ml to M3 the sup- 
pression is, again nearly exponential, but seems to cut off 
for M4. 

Due to the presence of the imposed azimuthal field, an 
additional anisotropy is induced to the system. The az- 
imuthal magnetic field is excerting a force opposing the 
movements in the directions perpendicular to it, which is 
seen as the suppression of the latitudinal and vertical ve- 
locity components, whilst the azimuthal velocity compo- 
nent is less affected (see Fig.ll6lpanel on the left). For the 
largest field strength (run M4) the azimuthal velocity ex- 
ceeds the other components, indicating that motions occur 
mostly along the field lines in contrast to the non-magnetic 
case (MO) in which the vertical (radial) direction was pre- 
ferred. This anisotropy is also seen in the quantities Ah 



and Av listed in columns 11 and 12, defined by Eqs. H47|l 
and 148() : the former is obtaining larger positive values as 
the magnetic field strength grows, the latter smaller neg- 
ative values and finally changing sign when the azimuthal 
motions grow larger than the radial ones. The fiuctuating 
magnetic field components (see Fig.ll6lpanel on the right) 
generated by the convective motions rather closely refiect 
this anisotropy bar the run M4, so that, e.g., the velocity 
component containing the most of energy generates the 
strongest fiuctuating magnetic field in its direction. 

The generated total turbulent stresses T^ , plotted in 
Fig. El show rather an unexpected behaviour. Although, 
in general, the fiuctuating velocity field becomes sup- 
pressed, the turbulent correlations are rather enhanced 
in magnitude when the strength of the imposed field is 
increased. In the hydrodynamic case MO, the stress com- 
ponents responsible for the generation of differential ro- 
tation, (Tr^) and (le^), are negative throughout the do- 
main, the former being the strongest in magnitude. For 
this component the profile retains its shape in the mag- 
netic regime, only the magnitude increasing as magnetic 
field grows. The growth is accompanied by the increase of 
the mean latitudinal flow (ux), so that it is roughly twice 
as strong in the case M4 than in MO. The two weak nega- 
tive peaks seen in Tg^ develop in a less systematic manner, 
the one at the interface of the convective and overshoot 
layers remains more or less unaffected by the growing mag- 
netic field, but the one near the surface changes its sign 
and gains in magnitude. The behaviour of T^g does not 
show any clear trend, but changes its magnitude, profile 
and sign rather irregularly. The mean azimuthal fiow, gen- 
erated by the rO component, also lacks a systematic be- 
haviour, although it tends to change its sign averaged over 
depth, being preferably positive in the regime of strong 
magnetic field. 

The coefficients -fftot, ^ot, and Mtot, calculated from 
the horizontal averages of the total stresses, are plotted 
in Fig. El as functions of the normalised field strength 
(3rn = \B\/ \/Jiopuf, where \B\ ~ Bo + bt- In the hydrody- 
namic case MO (see FigEI fourth row) the horizontal A- 
coefficient is positive, whilst the vertical coefficient and the 
quantity M are negative. As the magnetic field is applied 
to the system, the horizontal coefficient changes its sign 
and increases monotonically in magnitude (the dashed line 
in Fig. 1181) . the growth being somewhat reduced in the 
run M4. The vertical coefficient (the solid line) retains 
its negative sign, and also shows a monotonic increase in 
magnitude, with less evident signs of a suppression of the 
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Fig. 19. Volume averaged eddy heat flux tensor compo- 
nents Xrr/xt (left), and xor/xt (right) as functions of lat- 
itude for the Ra = 2.5 ■ 10^ runs. The errorbars denote 
the modified mean error of the mean as calculated from 
Eq. PHl. 



growth in the strong field limit. No obvious trend is vis- 
ible for the quantity M (dashed-dotted line), as already 
implied by the irregular behaviour of the stress (T^^). 
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Fig. 20. Same as Fig, [H but for tlie Ra = 10^ sets. 



As a comparison, we have overplotted tlie quenching 
functions for the weak field limit ^ 1 derived by 
Kitchatinov et al, H1994b|l in Fig, UHl (dotted lines). Their 
derivation is made in the limit of slow rotation for ar- 
bitrary field strength, due to which the results are not 
very well compatible with our moderate rotation case with 
Co = 2, but nevertheless provide useful comparison with 
the analytical turbulence models. In the hydrodynamical 
slow rotation limit, the vertical coefficient is predicted to 
be positive (see the dashed lines plotted in the first row 
of Fig, m or Fig, 2 of KR93), whilst the horizontal co- 
efficient vanishes. The inclusion of a weak magnetic field 
causes the appearance of an additional horizontal compo- 
nent negative in sign, which is enhanced as the magnetic 
field grows (oc 0^), except for the strong field regime, 
where the coefficient becomes quenched (oc Pm^)- In our 
case, the horizontal component is clearly enhanced, even 
so that the growth exceeds the prediction of the analyti- 
cal model (plotted on top with the dotted line, scaled to 
match the sign and hydrodynamical value of the calcu- 
lated i7-profile) and also shows a sign change compared 
to the hydrodynamic case, but the signs are reversed in 
comparison to the analytical results. Although the growth 
becomes somewhat suppressed in the strong field limit, it 
is significantly milder than predicted by the analytical in- 
vestigation. A quenching function of a similar form (plot- 
ted on top of the calculated F-profile with a dotted line, 
scaled to match the hydrodynamical value) was derived 
for the vertical component, accompanied by a sign change 
from positive to negative with increasing magnetic field. 
In our calculations the vertical component remains neg- 
ative and increases in magnitude as function of /3,„. The 
growth is also much stronger than the one predicted by 
the analytical expression and no significant quenching, as 
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expected from the analytical theory (oc /S^"^), however, is 
found in the strong field regime. 

In summary, the produced set of magnetic runs with 
an imposed azimuthal field indicate that the 9(j) and re/) 
components of the turbulent stresses, capable of gener- 
ating horizontal and vertical differential rotation, respec- 
tively, increase in magnitude with increasing imposed field 
strength whilst the rO component, capable of affecting 
the meridional circulation pattern, shows no clear trends 
although also gains its maximum value for the case of 
the strongest magnetic field. The model setup, however, 
hardly represents the situation in a real stellar convection 
zone, where magnetic fields arise due to the turbulent con- 
vective and large-scale shearing motions, exhibiting gra- 
dients and thereby mean electric currents leading to the 
Malkus-Proctor effect, which has been excluded from this 
investigation. 

5.8. Turbulent heat transport 

In Figs. El and [201 we show the latitudinal variation of 
the volume averaged eddy heat flux tensor components 
Xrr and xsr, defined via Eq. (021 ■ We use an estimate of 
the turbulent conductivity Xt = as normalisation. The 
error is estimated as in Figs. 1131 and 1141 

The radial component Xrr is always positive, indicat- 
ing that the heat flux is directed radially outwards, except 
at the equator for the case of highest rotation (ColO-00). 
In this calculation the temperature stratification becomes 
distorted by the strong rotational influence so that 5r be- 
comes negative. The latitudinal trends for clear: 
maximum occurs at the equator with diminishing magni- 
tude towards the pole and with increasing rotation rate. 
There seems to be, in general, a tendency for the mini- 
mum to occur at a latitude 9 « —45° ... — 60°, indicating 
a minimum in the efficiency of convection. These trends 
concur with the 'deeper layer' results of Pulkkinen et al. 
lfT(M|l . However, the occurence of the minimum in the 
turbulent heat flux at mid-latitudes seems to contradict 
the one obtained from the Nusselt numbers (see Fig. |HI), 
which indicates that the efficiency of convection decreases 
monotonically from the equator towards the pole. This 
discrepancy can be explained as follows: the value of the 
Nusselt number, Eq. (|^ . depends on only f'^'^^ = K^f^, 
where the difference A is taken using the values at the 
layers zi and Z2. In other words, the temperature strati- 
fication is assumed to be linear. If this was true, the ra- 
dial superadiabacity, (5^, would be proportional to 
and the latitudinal profiles of Xrr would show the same 
qualitative trend as the Nusselt number since {u'^T') itself 
and Nu show similar tendencies as function of latitude. 
The superadiabatic layer, however, gets shallower from 
the equator to the pole, which affects immediately 
because it depends only on the values of e at depths zi 
and Z2. On the other hand, the convection zone is more 
superadiabatic (i.e 5r is larger) at mid-latitudes than near 
the poles. Since Xrr takes into account the whole strati- 



fication rather than just the endpoints, it probably gives 
a more reliable estimate of the convection efficiency than 
the Nusselt number does. 

The xor component shows a less coherent picture. For 
slow and moderate rotation xsr is negative throughout 
the latitude range, the maximum being at low latitudes 
(6 ~ —10 . . . — 30°). For the case Co4 xer obtains positive 
values near the pole and for the higher rotation cases there 
is a definite positive peak at high latitudes, and another 
near the equator. 

Due to the boundary conditions we cannot actually 
have a temperature difference between the top layers of 
the boxes situated at different latitudes, but the radial 
gradient and the value at the bottom of the box are not 
constrained. In accordance with the latitudinal distribu- 
tion of Xrr, we find that the temperature gradient in the 
equatorial cases is always shallower than at higher lati- 
tudes, and that the mid-latitudes do seem to be the least 
favoured in the sense that there the gradient is the steep- 
est. If this situation was naively moved onto a spherical 
star where we assume the core to have a constant surface 
temperature below the convection zone, this would result 
in a warm equator and pole, and a cooler region at mid- 
latitudes. However, this approach is probably too simple 
since it does not take into account the latitudinal trans- 
port which shows a more complex pattern. For slow rota- 
tion (CoOl to Col), the negative values of XrS indicate an 
equatorward heat flux at all latitudes. In the intermediate 
regime (Co2 and Co4) there is a poleward flux at the polar 
regions and a less pronounced equatorward one near the 
equator, supporting the view of the cooler mid-latitudes 
discussed above. This tendency, however, is reversed in the 
rapid rotation regime, so that a poleward flux appears at 
the equatorial regions, and an equatorial one at high lati- 
tudes, which kind of a flux would act as smoothing out the 
latitudinal temperature difference generated by the radial 
heat flux. 

6. Conclusions 

We have performed local three-dimensional hydrodynam- 
ical calculations of compressible convection with the in- 
clusion of an overshoot layer in Cartesian volumes placed 
at various latitudes and with largely varying rotation. We 
find that the convectivc structure changes from an ordered 
cellular pattern of slow and moderate rotation (Co < 2) 
to a highly irregular small-scale pattern for rapid rota- 
tion (Co > 4). The efficiency of convection is noted to, in 
general, decrease as function of increasing rotation, while 
the latitudinal dependence indicates a minimum at mid- 
latitudes and more efficent convection at the equator and 
near the poles. 

The kinetic helicity, 7i, shows a similar trend as 
previous numerical studies of magnetoconvection (e.g. 
Ossendrijvcr et al. ,2002) and helioseismic inversions 
(Duvall & Gizon l^nUnjl or large positive values at the south 
pole and an approximately cos 0- latitude dependence for 
slow and moderate rotation. We, however, find that for 
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rapid enough (Co > 4) rotation the sign of the hehcity 
changes and the latitudinal distribution follows more a 
sin ^-dependence at latitudes Q < —15° than the cos9- 
profile typical for the two aforementioned studies. It has 
been concluded that the a computed from magnetocon- 
vection calculations qualitatively agrees with the first or- 
der smoothing result a « —^t'TI (e.g. Ossendrijver et 
al. I2001|l . If this behaviour carries over to the rapid rota- 
tion regime, the a-effect would also be subject to a sign 
change and simultaneous decrease of magnitude at higher 
latitudes. 

The convective overshooting is reduced monotonically 
as a function of rotation. The latitudinal dependence 
shows two regimes; (i) for slow and intermediate rotation 
the overshooting depth decreases as latitude increases, 
and that (ii) for rapid rotation the trend is the oppo- 
site. Imposing an azimuthal magnetic field causes the over- 
shooting depth to decrease monotonically as a function of 
increasing field strength. 

The mean flows, generated by the Reynolds stresses, 
also show two distinct regimes; in cases of slow and in- 
termediate rotation the latitudinal flow is poleward and 
the azimuthal one retrograde. Both components exhibit a 
more or less linear shear profile in the convectively unsta- 
ble region. For rapid rotation, the signs of the mean flows 
change so that there is a strong equatorward flow in a 
narrow region near the surface and a weak poleward flow 
in the bulk of the convection zone. The azimuthal flow is 
prograde in the bulk and with a retrograde dip near the 
upper boundary. 

The volume-averaged Reynolds stress components, 
{Qe4>)v and {Qr4>)v, responsible for the generation of hor- 
izontal and vertical differential rotation, respectively, indi- 
cate that the angular momentum is transported equator- 
wards in all cases where rotation is significant (a positive 
horizontal A-effect). The radial transport occurs inwards 
(negative vertical A-effect) for slow and moderate rota- 
tion and towards the surface for rapid rotation (positive 
vertical A-effect). The latitudinal distribution shows that 
the inward transport in situated mostly at high latitudes 
and the outward transport near the equator. The merid- 
ional Reynolds stress component, {Qre)v is negative for 
slow rotation and reaches positive values at high latitudes 
near the pole. We find that both A-effects are subject to 
rotational quenching, indicating the suppression of differ- 
ential rotation in rapidly rotating stars. Comparison to 
the analytical turbulence model of KR93, based on the 
quasi-linear approach, shows a rather close agreement for 
the horizontal A-effect near the polar regions, whilst the 
values obtained from our calculations are systematically 
larger near the equator and less suppressed as the rota- 
tion is increased. The regime of positive vertical A-effect 
for slow rotation is not found from our calculations, but for 
intermediate rotation the agreement with KR93 is fairly 
good. For the rapid rotation regime, however, the positive 
vertical A-effect obtained in this study is in disagreement 
with the analytical theory. 



We have also produced a set of magnetohydrodynam- 
ical calculations by imposing a mean azimuthal magnetic 
field on a saturated convection pattern with Co = 2. As 
the imposed field does not have any gradients, no mean 
currents are present in the system, due to which the large- 
scale Lorentz force vanishes. The action of the convective 
motions on the mean azimuthal field, however, leads to 
the generation of a fiuctuating component of the mag- 
netic field and non-zero Maxwell stresses. In the magneto- 
hydrodynamical regime both the Reynolds and Maxwell 
stresses contribute to the angular momentum transport; 
the magnitude of the total turbulent stresses is observed 
to increase monotonically, although the turbulent velocity 
field, in general, becomes suppressed. The horizontal A- 
effect is observed to change sign from positive to negative 
as the field strength increases, indicating a new regime of 
poleward transport due to the presence of the magnetic 
field. 

The turbulent heat transport indicates a rather strong 
latitudinal dependence of the radial heat flux, whereas 
the latitudinal flux shows no distinct tendencies. A simple 
analysis of the radial heat fiux indicates that the equato- 
rial regions should be the warmest, followed by the polar 
regions, and a cool belt at mid-latitudes would appear. 
However, this interpretation is most probably based on 
too simple assumptions since the thermal structure is also 
affected by the latitudinal heat fiux and the meridional 
circulation. On the other hand, the latter depends on the 
thermal structure so a straightforward analysis as above 
is inadequate. 

Having the knowledge of the Reynolds stresses and tur- 
bulent heat transport as functions of rotation and latitude 
it is, in principle, possible to investigate the implications of 
these local results on the global dynamics (rotation law, 
thermal structure) by utilizing a mean-field model in a 
spherical geometry. However, this analysis is not within 
the scope of this paper. 
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Fig. 12. A polar representation of the horizontally averaged Reynolds stress components {Q04,) (top), {Qr4,) (middle), 
and {Qro) (bottom), normalised by the turbulent viscosity. Data from the sets Col (left) and ColO (right) is shown. The 
dashed lines denote the top and bottom of the convectively unstable layer, and the dotted lines indicate the latitudes 
for which data exists; to represent the latitudinal dependence of the stresses smoothly, a cubic spline interpolation 
between the datasets has been made. Note that the stresses are presented as they would appear in spherical coordinates 
on the southern hemisphere. 



